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Abstract 


The  aim  of  this  paper  is  to  study  the  convergence  properties  of  a  Time 
Marching  Algorithm  solving  Advection-Diffusion  problems  on  two  domains 
using  incompatible  discretizations.  The  basic  algorithm  is  first  presented, 
and  theoretical  or  numerical  results  illustrate  its  convergence  properties. 
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1  Introduction 


Domain  decomposition  methods  have  become  an  efficient  strategy  for  solv¬ 
ing  large  scale  problems  on  parallel  computers  ([1],  [2],  [3],  [4],  [5],  [6]). 
Nevertheless,  they  can  also  be  used  in  order  to  couple  different  models  [11], 
[18],  [19]  and  [21].  In  this  paper  we  will  examine  a  domain  decomposition 
strategy  which  can  be  applied  to  such  situations. 

This  approach  was  introduced  in  order  to  solve  several  difficulties  that 
occur  in  fluid  mechanics.  In  particular,  our  aim  is  to  introduce  several 
subdomains  in  order  to  locally  introduce  an  enriched  model  next  to  a  domain 
boundary.  For  this  purpose,  we  propose  to  fully  overlap  the  subdomains 
and  to  couple  the  solutions  through  natural  “friction”  (Neumann)  forces 
acting  on  the  internal  boundary  of  the  domain,  these  friction  forces  being 
updated  inside  the  time  marching  algorithm  used  for  the  solution  of  the 
initial  problem. 

The  theoretical  study  of  our  method  will  be  done  on  an  Advection- 
Diffusion  problem,  which  will  serve  as  our  model  problem  from  now  on. 
The  analysis  will  be  made  at  the  continuous  level,  independently  of  any 
discretization  strategy,  which  means  that  the  derived  results  will  be  mesh 
independent.  The  use  of  friction  (Neumann)  coupling  boundary  conditions 
makes  the  convergence  analysis  somewhat  different  of  the  analysis  done  in 
Kuznetsov  [22]  or  Rannacher  [23]  in  their  study  of  explicit  Schwarz  additive 
methods  for  time  evolution  parabolic  problems. 

In  the  next  section  we  will  describe  this  model  problem.  In  the  third 
section  we  will  present  our  algorithm  for  some  basic  cases.  The  fourth  section 
will  treat  the  one-dimensional  stationary  problem.  We  will  show  also  that 
the  convergence  of  this  method  can  be  improved  by  introducing  a  relaxation 
parameter  [5].  The  fifth  section  will  focus  on  the  linear  convergence  of  the 
implicit  version  of  the  coupling  algorithm  in  the  general  multidimensional 
case.  In  the  last  section  we  study  the  numerical  stability  of  the  explicit 
algorithm.  Practical  applications  of  the  proposed  algorithm  to  real  life  CFD 
problems  can  be  found  in  [14],  [19],  [20],  and  [21]. 

2  The  Model  Problem 

Consider  a  bounded  domain,  ft  of  Rn  such  that  its  boundary  <9ft  is  lips- 
chitzian  and  ft/oc  a  connected  domain  of  Rn  with  ft;oc  C  ft  (fig.  1).  The 
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Figure  1:  Description  of  the  computational  domain. 

boundaries  of  the  two  subdomains  are  defined  as  follows: 

rfc  =  dQ,  D  dftioc,  ( internal  boundary) 

Ti  =  dttioc  n  17,  (  interface) 

Toe,  =  dtt\Tb.  (farfield  boundary) 

We  denote  by  n  the  external  unit  normal  vector  to  dQ,  or  di 1ioc- 
We  will  make  use  of  the  following  notation 


IMIo.o  =  IMIl2(o) 

IbfUo  =  Ibll^(o) 

Mi,o  =  ||Vv||z,2(o) 

where  0  is  an  open  bounded  domain  of  Rn. 

Let  v  be  the  velocity  field  inside  a  given  incompressible  flow  such  that: 


divt;  —  0  in  0, 
v.n  —  0  on  IV 


(i) 


We  consider  the  following  convection- diffusion  model  problem: 
Find  '~p.  a  real  valued  function,  defined  on  J7  and  satisfying 
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div(v<p)  —  uA(p  ~  0  in  fi, 

<  <t>  =  <f°°  on  Too,  (2) 

<p  =  0  on  rv 

Above  v  is  the  flow  velocity  and  v  is  the  diffusion  coefficient.  Problems  of  this 
form  typically  occur  in  fluid  mechanics,  gas  dynamics  or  wave  propagation. 

Most  CFD  algorithms  will  in  fact  consider  the  solution  of  this  problem 
as  the  stationary  solution  of  the  evolution  problem  (3)  described  below  : 


Find  6  :  ft  x  (0,  T)  — *■  R  such  that, 


^  +  div(i><p)  -  vA 4> 


0  in  Cl  x  (0,T), 


<j)  =  r  on  Tco  x  (0,T), 

<  (3) 

4>  =  0  on  Tji  x  (0,T), 

„  </>(0)  =  <f>o  in  ft. 

The  general  CFD  algorithm  consists  then  in  integrating  (3)  in  time  until 
reaching  a  stationary  solution. 


3  General  Algorithm 

3.1  Time  Continuous  Case 

Let  us  introduce  the  local  subdomain  ft;oc  (see  fig.  1)  which  has  as  external 
boundary  and  let  us  consider  the  trace  <t>ioc  of  <j>  on  the  subdomain  ft;oc, 
as  an  independent  variable,  to  which  we  associate  an  arbitrary  independent 
initial  value  <f>0i  ^  4>o\uloc-  We  now  replace  the  evolution  problem  (3)  by  the 
following  evolution  system  : 

Find  <f>  (resp.  <j>ioc)  :  ft  -*■  R  (resp.  ft;oc  —  R)  satisfying 

+  div(u^)  —  uA<j>  =  0  in  ft  x  (0,T), 

<  <t>  =  <t>°°  on  rTO  x  (0,T),  (4) 

d<f>  df>ioc  . 

1  =  "“^T  on  riX(0’r)' 
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(5) 


+  di v(v<f>loc)  ~  vA<f>ioc  =  0  in  ft/oc  x  (0,T), 
at 

<  <t>loc  -  0  on  r6  x  (0,  T), 

<f>ioc  =  <f>  on  Ti  x  (0,  T), 

(f)(0)  =  4>o  in  ft,  <f>ioc( 0)  =  <f>oi  in  loc ■  (6) 

Remark  3.1  The  global  problem  (4)  with  the  initial  condition  (6)  has  no 
no-slip  boundary  condition.  This  suppresses  the  boundary  layer  which  ap¬ 
pears  at  low  viscosity  and  facilitates  the  numerical  solution  of  this  problem. 
The  boundary  layers  are  modeled  by  the  local  problems  (5) -(6)  which  are 
only  to  be  solved  on  a  small  domain  f 1loc,  with  a  very  fine  discretisation  if 
needed.  The  two  problems  are  only  coupled  by  their  boundary  conditions  and 
not  by  volumic  interpolation. 


3.2  Time  Discrete  Case 

The  general  algorithm  that  we  propose  for  the  solution  of  our  model  problem 
(2)  is  as  usual  to  integrate  in  time  the  evolution  problem  (4)-(5)-(6)  until  we 
reach  a  stationary  solution.  This  integration  in  time  is  then  achieved  by  the 
following  uncoupled  semi-explicit  algorithm,  where  the  operators  are  treated 
implicitly  inside  each  subdomain  and  where  one  of  the  coupling  boundary 
conditions  is  treated  explicitly  and  the  other  is  treated  implicitly: 

•  set  4>°ioc  -  <t>0i  and  <t>°  -  <t>o, 

•  then,  for  n  >  o  A  and  4>n  being  known,  solve  successively 


r  ct1  - 1 


At 


—  +  divO^1)  -  vAff*1 


ct1 

ct1 


0  in  fl/oci 
<f)n  on  r,-, 
0  on  Ti, 


(ffi+l  _ 

At 


+  div(u^n+1)  -  vA <f>‘ 


kn+ 1  — 


0  in  ft, 


<f>n+1 

d(f>n+1 

dn 


4>°°  on  Too, 
—  y - -M—  Dn 


dn 


(7) 


(8) 
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Remark  3.2  We  have  a  full  decoupling  between  (7)  and  (8).  They  can 
(and  actually  will)  be  discretized  and  solved  by  two  independent  solution 
techniques . 

Remark  3.3  The  fully  implicit  version  of  this  method  consists  in  replacing 
the  condition  : 

cf  =  <t>n  on  rt 

by  the  condition  : 

4>lV  =  ^"+1  on  T{. 

The  two  subproblems  are  then  coupled  at  each  time  step . 

Remark  3.4  If  we  replace  in  (8)  SI  by  SI e  defined  as  follows  : 

SIe  —  Sl\Slioc , 

and  Tb  by  and  if  we  set  At  =  oo,  we  obtain  a  nonoverlapping  version 
of  our  strategy ,  which  is  a  standard  Dirichlet- Neumann  algorithm  [16],  [17] 
and  therefore  requires  a  relaxation  strategy  to  converge . 

Remark  3.5  The  initial  condition  <f>0i  is  not  assumed  to  be  equal  to  <f> o  on 
the  local  subdomain  Slioc  because  in  most  cases  this  condition  is  impossible 
to  impose  at  the  discrete  level  since  the  grid  used  on  Slioc  will  be  in  general 
different  from  the  grid  used  on  SI.  In  addition ,  even  if  we  assume  <f)0t  =  <j>o, 
we  will  not  have  <f>foc  =  <j)n  on  Slioc  unless  we  use  the  fully  implicit  algorithm 
on  compatible  grids . 

4  Stationary  one-dimensional  case 

For  At  =  +oo,  the  above  algorithm  can  be  written  : 

•  set  (f>°oc  =  <f>o  and  <j>°  =  <f>0, 

•  then,  for  n  >  0,  <f)foc  and  d>n  being  known,  solve 

div(u<^J^1)  —  vA4>™£}  =  0  in  fyoc, 

C^1  =  <t>n  on  Ti,  (9) 

,  Ct1  =  0  on  Tb, 
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=  0  on  ]/i2,l[, 

9f\k2)  =  9!"-1)(A2),  (13) 

*4”’(i)  =  »> 

=  0  on  ]0, 1[, 

9«”»(0)  =  a,  (14) 

*4n)'(  1)  =  4n)'(i)- 

By  introducing  two  relaxation  parameters  and  62,  we  can  also  intro¬ 
duce  the  following  variant  of  the  above  algorithm  : 
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(15) 


(")'  («)" 
V<p2  -V(p\> 

» 


0  on  }h 2, 1[, 


(1)  =  b, 

2)  =  ^2^1"  ^(^2)  +  (1  —  ^2)^2  ^2)? 

(nY  (n)n  n  In  1  r 

v<p\  —  v<p\  =  0  on  JO,  1[, 

rfn)(0)  =  a,  (16) 

^n)'(l)  =  + 

We  shall  now  exhibit  the  conditions  under  which  the  algorithm  (15)-(16) 
converges,  and  those  for  which  this  convergence  is  optimal.  For  this  purpose, 
we  write  the  interface  solution  under  the  form 


*,<">'(1)  =  *>'(!)  + 7".  (1?) 

4“’(fc2)  =  V(l>2  )+<”.  (18) 

where  is  the  solution  of  the  initial  problem  (12).  Using  the  analytical 
solutions  of  the  problems  (15)  and  (16),  we  obtain  the  following  induction 
formula 


with 


=  Min 


\ 

-y(Tl-l)  )  ' 


(19) 


(  1-02 


Min  = 


1  *i(l  -02)0 
\  -  1 


eS-e-d)^  _  1) 

V 


-  1) 

(e(p)(^2-i)  _ 


+  (l~0i) 


(20) 


This  iterative  process  converges  if  the  spectral  radius  of  the  matrix  Min  is 
less  than  1.  A  direct  but  tedious  calculation  then  yields: 


Lemma  4.1  The  spectral  radius  of  the  transfer  matrix  of  the  algorithm 
(15)-(16)  is: 
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(21) 


p(M1N )  =  max[^\D  ±  \[D2  -  AR\] 

with 

D  =  2  -  (»i  +  <h)  +  Z  i  (22> 

i?  =  (l-^)(l-^)-  (23) 

From  this  calculation  we  obtain  the  following  results: 

i)  When  /i2  goes  to  1  (nonoverlapping),  D  goes  to  +oo,  and  then,  p(Min) 

goes  to  +oo.  There  is  no-convergence  at  this  limit. 

ii)  The  optimal  convergence  is  obtained  in  the  case  where  all  the  eigenvalues 

of  the  matrix  Min  are  zero,  i.e.,  when  :  D  =  0  and  R  =  0.  The  latter 
conditions  imply  in  particular 


0\  =  1  or  $2  =  1. 

If  we  choose,  in  addition,  0\  -02,  the  condition  D  -  0  implies  h2  =  0. 
In  this  case  the  subdomain  Q/oc  is  equal  to  the  whole  domain,  and  the 
associated  algorithm  is  no-longer  of  interest. 

iii)  The  convergence  of  the  method  depends  symmetrically  on  both  relax¬ 
ation  parameters. 

According  to  ii)  it  is  reasonable  to  take  one  of  the  0{  equal  to  1  and  call  the 
other  0. 

By  setting: 

^  ~  1  e(v/i/)(h2-l)  _  i 

we  then  have 

p(MIN)  =  \1-0A\. 

In  this  case,  setting 

a  _  /i  _  e(-^)(e^)fe2  -  1)  x 
v0pt  -  IT  e(v/v)(h2- 1)  _  l  J  ’ 

which  is  <  1,  we  get  the  following  convergence  results: 


(25) 

(26) 


,(v/v)h.2 


1) 


(<1A\ 
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Theorem  4.1  1)  The  convergence  is  optimal  (convergence  in  1  iteration) 

if 

e  =  6opt.  (27) 

2)  The  algorithm  converges  for  all  6  in  ]0.  -j  [. 

Corollary  4.1  1)  The  case  without  relaxation  (6  =  1)  converges  only  if : 


i.e.,  by  setting  d  =  1  -  h2  (overlapping  length),  only  if  : 

d  >  -Log- - - — j—  (stability  condition). 

v  ( 1  +  e~v lv) 

2)  When  v  goes  to  zero,  we  must  have  d  > 

Remark  4.1  This  theorem  states  that  the  application  of  the  algorithm  (15)- 
(16)  to  the  time-independent  problem  (12)  converges  only  if  the  overlapping 
d  is  sufficiently  large.  In  the  same  situation,  we  will  see  that  if  the  problem 
(12)  can  be  regarded  as  the  steady  solution  of  a  time-dependent  problem 
and  we  apply  our  strategy  to  this  evolution  problem,  the  resulting  algorithm 
will  converge  to  the  same  steady  solution  but  with  less  restrictions  on  d. 
This  motivates  the  introduction  of  the  time  marching  algorithm  of  section  3. 
Moreover,  this  time  marching  technique  is  well  adapted  to  nonlinear  problems 
such  as  those  encountered  in  fluid  mechanics  (see  [If],  [19],  [20],  and  [21]). 


5  Implicit  Time  Discretization 

5.1  The  General  Algorithm 

This  section  deals  with  the  convergence  analysis  of  the  proposed  algorithm  in 
multiple  dimensions  when  one  uses  the  fully  implicit  version  of  our  strategy 
(4)-(6)  : 

•  Set  <b°loc  =  <t>0i  and  <j>°  =  <p0; 

•  then,  for  n  >  0 ,<j>?oc  and  on  being  known,  solve 
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<j>- n+1  -  (f)71 


+  div(v<f)n+1)  —  vA<f>n+l  =  0  in  ft, 


dcj>n+1 


—  4 oo  on  r  oo  j 

=  Jf-  o„  r(, 

on 


f  i7l-\-l  _  in 

'K  At  +  div(<,^‘)  -  <'M£l  =  0  “  «!«=> 

<  f;*'  =  r*1  on  i„  <29) 

Ct1  =  0  on  r4. 

5.2  Convergence  Analysis 

Before  establishing  the  convergence  result  we  shall  state  the  preliminary 
results  that  are  central  to  the  proof  of  the  convergence  of  our  algorithm. 
The  first  result  states  the  basic  L2  and  H1  local  estimates. 

Lemma  5.1  We  have  the  following  estimates: 

^nv+l  -ct'iioa.  +  4rK-«t1ii,sw 

<  -  CJIww  (30) 


Un+1-&lo,x  <  1  +  2,Ate"^-^llW.  (31) 


ir+1-OS, 


<w  £  (TT2^)”+,|I/-&|I (32) 


l®"+1  -  cyilo.o,,.  +  2nAfX>i+1 


<  11^  -  tfjlww  <  », 

where  c  is  the  Poincare  constant  on  subdomain  tlioc- 
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Proof  of  lemma  5.1 


Substracting  (29)  from  (28),  multiplying  the  result  by  <f>n+1  -  Ct*  and 
integrating  by  parts  over  f l;oc,  we  obtain  the  classical  following  relation: 


By  using  the  Cauchy-Schwarz  inequality,  we  obtain  the  estimate  (30).  The 
second  estimate  (31)  follows  by  using  the  Poincare  inequality  with  c  the 
Poincare  constant  bounding  the  squared  H 1  seminorm  of  any  function  v  of 
•ffH^/oc)  with  zero  trace  on  I\-  by  its  squared  L2  norm.  By  induction  we 
also  obtain  the  basic  L2  estimate  (32).  And  finally,  we  obtain  the  estimate 
(33)  by  summing  (30). 

■ 

The  above  lemma  states  that  the  restriction  of  <f>n+1  -  Ct*  to  f 1ioc 
converges  to  0  in  both  L 2  and  H1  norms.  We  shall  establish  now  other  L 2 
and  H 1  local  estimates.  Let  6xn  be  defined  by 


8xn  = 

and  let  G  be  defined  by 


(^n+1-^t1)-(^n-Cc) 

At 


G(n ) 


Un  -  Cell 


2 

0 


oc 


+  IC-Cd 


2 

1 J ^loc  * 


(35) 


(36) 


Lemma  5.2  We  have  the  following  estimates: 


<  |j(G(n)-G(n  +  l))  (37) 

G(n  +  1)  5  +  (38) 

Proof  of  lemma  5.2 

Substracting  the  two  first  equations  in  (28)  and  (29),  multiplying  the  result 
by  8xn  and  integrating  over  Qioc  we  obtain 
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/ 

\6x 

n|2+  / 

div(„(^n+1  -  4>^cl))Sxn 

j^loc 

JO'loc 

+v 1 

1 

V(d>n+1  - 

J 

&loc 

-VJ 

f 

dQloc  on 

Using  the  boundary  conditions  in  (28)  and  (29),  and  the  Cauchy-Schwarz 
inequality  we  obtain 


IMS*,.,  <  jii>Cr+1  -  OW  +  jMW 

+^t\r  -  -  ^”+1  -  <C'\G 


(39) 


2  At' 


Using  now  the  relation  (30)  (lemma  5.1)  leads  to  the  first  estimate  of  our 
lemma.  In  fact,  this  estimate  implies  that  G  is  a  decreasing  function.  This 
property  then  yields 


n+1 


(n  +  2  -  p)G(n  +  1)  <  £<?(*) 


i—p 

n+l 


iii-ii?  "+i 


<  X)  “  $ocllAoc  +  11^*  “  ^/ocIloA 


i=p 


t=p 


(40) 

Using  again  the  relations  (30)  and  (31)  (lemma  5.1)  yields  the  second 
estimate  (38).  And  the  lemma  is  proved. 


We  shall  establish  now  the  global  L2  and  H1  estimates.  Let  <£n+1  be 
defined  as  follows 


4>n+1  = 


<f>n+1  -  <p  in 

Ct1  ~  V9  in  ttioc, 


(41) 


with  the  solution  of  the  stationary  problem  (2).  By  construction, 
satisfies  the  following  equations: 
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f  ^+i-i^  +  div(^+1)-^+1  =  0  in  Shoe  U  (ft\n/oc), 

4>n+1  =  0  on  dn, 

4>n+1  continuous  across  T;. 

Let  A,  B\ ,  and  B2  be  defined  by  the  following  relations 

1 


A  = 


1  +  ( vc  —  2c\)At  ’ 
Bi  =  ^-A, 


(42) 


Hull2 

B2  =  At(]±-fZ  +  v)A, 

ci 

where  c  is  the  Poincare  constant  and  c\  >  0  is  an  arbitrary  constant  as  will 
be  seen  in  the  proof  of  the  following  lemma. 


Lemma  5.3  li  e  have  the  following  estimates 

1 

K2At  "u’“  ‘  2,r  11,14  “  2At 


(43) 


ll#n+,llS,o  <  A||#”||S, a  +  Bil|«*lSllw  +  B 2|<A”« 


uloc 


-r+1\l*„ 


(44) 


Proof  of  lemma  5.3 


Multiplying  the  equation  (42)  by  4>n+1  and  integrating  by  parts  over  ft/oc 
and  n\fi/oc,  and  taking  into  account  the  boundary  conditions  in  (42)  we 
obtain  the  following  relation: 


/  + v  [  \v4>n+i\2 

Ju  A  t  Jq 

/n  -  r+1¥“+1  =  0. 


(45) 
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On  ft/oc,  d>^cl  -  <f>n+1  satisfies  the  following  equation 

(Ct1  -  ^+1)  -  Mk  ~  fl  +  div[u(^+1  -  d>n+1)\  -  -  <T+1)  =  .0. 

_  (46) 

Therefore,  multiplying  the  above  equation  by  <fn+1 ,  integrating  by  parts  and 
using  the  Cauchy-Schwarz  inequality  we  obtain 


with  ci  >  0  arbitrary.  Combining  the  above  inequality  with  (45),  bounding 
the  local  norm  \f\i,Qloc  by  \f\t,n  and  using  the  Cauchy-Schwarz  inequality  we 
obtain  the  estimate  (43).  The  estimate  (44)  results  immediately  from  the 
estimate  (43)  by  applying  the  Poincare  inequality  on  Cl  with  c  the  Poincare 
constant.  And  the  lemma  is  proved. 

■ 

Finally,  we  are  in  a  position  to  state  the  main  result  of  this  section. 

Theorem  5.1  The  solution  of  the  algorithm  (28)-(29)  converges  linearly  in 
Hl( 0)  to  the  solution  of  the  stationary  problem  (2),  for  all  values  of  At  and 
all  choices  of  Choc. 

Proof  of  theorem  5.1 

Let  ci  be  chosen  such  that  vc  -  2 cj  >  0.  Using  the  relation  (44)  (lemma 
5.3)  we  obtain  by  induction 

-  r+l-‘\l,ai  j. 

Since  A  <  1  by  assumption  on  ci,  this  implies 
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lir+1|lo,Q  <  + Y,  Iloilo, !!,«  + 

i=n+l~p 

B 2  Y  WJ1  -  #'+1l?A«)- 

t=n4*l—  V 

Now,  using  (37)  (lemma  5.2)  and  (33)  (lemma  5.1)  we  obtain 


(48) 


Il^”+1|l5,a  <  A”||^+1-"||Sin  +  A(ft^(G(n+l-p)-G(n+l))+ 

ft  A 7ntf-'  -  r+i-i§,a,j. 


'2i/A*' 

The  same  relation  written  between  0  and  n  +  l—p  yields 


(49) 


||#"+1-’’IIS.0  <  A^UX^  +  ABi’-iGW-Gin-p  +  l)) 

+ft  2^Ai^°  ~  ^'“11°,""*)' 

By  combining  this  relation  with  (49),  we  finally  obtain 


IW”+1|lo,a  <  .4r-1|!o°llf1;.i> 

+  ^G(O)  +  ft  “  CIISaJ 

+  A(B^G(  n  +  1  -  p)  +  Ai11^1"’  “  ^+1_’II8a~)- 

Choosing  p  such  that  n  =  2p  +  q,  q  >  1  and  using  (38)  (lemma  5.2)  we 
conclude  that 

\\r+%n  <  An+1C2  +  AP+'Cs  +  CiWtioc  ~  (50) 

which,  from  (32)  (lemma  5.1),  implies  the  linear  convergence  of  ||^n+1||o,n 
to  0. 

On  the  other  hand  by  combining  (37)  (lemma  5.2)  and  (43)  (lemma  5.3) 
we  obtain 
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(51) 


1  +  §0> 


1 


n.  i|2 

1,12  ^  gAfllV  llo, 12 


'n+l12  -  <  t4- 


+  2^At{G{n)  -  G{U  +  1)}  +  (^f  +  2 : Moc  '  ^ 

Therefore  by  using  (30)  we  obtain 


(j^-oi)il^+,HS.n  +  |wn+1li<i5 
+S|A((G(“)"'S("  +  1))+(J2T 


V.Mfoc-  ^llo, 0(00  5 
+  2A  2vAt  J 


(52) 


Our  result  follows  then  from  (38)  (lemma  5.2),  (32)  (lemma  5.1),  and  the 
linear  convergence  of  ||<£"||o,fi- 


5.3  Convergence  of  a  Fixed  Point  Method  for  the  Implicit 
Scheme 

The  implicit  scheme  proposed  in  this  section  couples  the  global  and  the  local 
problem.  To  uncouple  them,  it  is  advisable  to  use  the  fixed  point  algorithm 
below  : 

•  set  </>?OC;0  =  iboi  and  (f>°  =  ^o, 

•  then  for  A;  >  0, 4>nk^  being  known, 
solve 


r  «S+i  - « 


loc 


At 


+  div(u<%^+1)  -  vA 4>^ocM\ 

=  0  in 

^"otfc+l 

=  rk+1 

on  r,-, 

=  0  on 

n, 

(53) 
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r  Cm  -  r 

At 


+  div(v^+j)  -  vA4>1+\  =  0  in  Q, 

<t>nkl\  -  0n  Too, 
vd&kW  / On  -  vd(t>iocMirdn  on  rf>- 


We  will  study  now  the  algorithm  (53)-(54).  By  setting 

i>loc,k,q  =  <C,fc+l  —  ^loc,q+ 1’ 

i>K,  =  (rk+'-r,+l\ 

we  observe  that  i>ioc,k,g  an^  i’k.q  verify  the  following  equations  : 


(54) 

(55) 

(56) 


("  Tploc,k,ql div(yiploc,k,q)  —  0  in  Tl\oc, 

1>locJt,q  =  Vfc-1,<?-1  On  I\,  (57) 

V* loc,k,q  =  0  On  IT*  f, , 


1pk,q/At  +  div(vipk,q )  -  vAi)k,q 

i>k,q 

d^k,q 

V - — 

dn 


=  v 


0  in  £1, 

0  on  Too, 

dibloc, k,q 


dn 


on  IV 


(58) 


If  At  is  sufficiently  small,  we  prove  in  [15]  that  ibk,q  and  ipioc,k,q  converge 
linearly  to  zero.  Hence  the  sequences  <f>^+1  and  are  Cauchy  sequences 
which  converge  linearly  towards  the  unique  solutions  <j)n+1  and  of  the 
implicit  scheme.  This  guarantees  the  convergence  of  the  above  fixed  point 
algorithm. 


6  Numerical  Analysis  of  the  Stability  of  the  Al¬ 
gorithm  (7)-(8) 

In  this  section  we  focus  on  the  application  of  the  explicit  time  marching 
algorithm  (7)- (8)  studied  in  the  previous  sections  to  the  numerical  solution 
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of  the  steady  problem  (2). 
r&  in  (8)  is  explicit 


We  first  assume  that  the  boundary  condition  on 


dn 


=  v 


dn 


) 


so  that  the  resulting  algorithm  is  parallel  (Jacobi  type). 

Here,  Q,  denotes  the  domain  surrounding  the  obstacle  (an  ellipse  in  our 
numerical  study)  as  described  in  Figure  1.  The  global  and  local  domains  are 
discretized  by  fully  overlapping  compatible  finite  element  grids.  The  global 
mesh  contains  1378  nodes  and  2662  elements  (see  figure  7).  Further  the  time 
marching  algorithm  is  being  initialized  by  setting  (f> o  to  zero. 

In  a  first  step,  the  velocity  field  is  obtained  by  solving  the  following 
inviscid  incompressible  flow  problem: 


divu  =  0, 

curl  v  =  0 

Voo  =  (1,0), 

v.n  =  0  on  the  body  Tj, 

with  a  first  order  finite  element  method  using  the  same  global  mesh. 

If  we  set  v  =  0,  the  algorithm  may  or  may  not  converge  depending  on 
the  values  of  vAt.  More  precisely,  we  observe  that  the  algorithm  converges 
linearly  when  vAt  <  a0  and  is  divergent  otherwise.  This  is  graphically 
shown  in  the  figures  (4-7)  where  the  values  of  are  plotted  versus 

the  iteration  count  n  for  vA t  equal  respectively  to  10-6,  10-1,  1,  and  10. 
Further,  when  the  velocity  is  taken  sufficiently  large,  the  algorithm  becomes 
unconditionally  stable.  In  particular,  the  initialization  of  our  algorithm  by 
<f>o  =  0  with  1 1  Vqq 1 1  =  1,  v  =  0.1  and  At  =  100  leads  to  a  converging  algorithm 

(fig-  8)‘  .... 

By  intuition  such  a  behavior  seems  natural.  An  overestimation  of  the 
solution  <j>n  at  the  interface  F,  implies  an  overestimation  of  the  friction  forces 
on  Tb.  For  sufficiently  small  time  steps,  this  overestimation  will  not  affect 
the  value  of  4>n+1  on  T;  and  can  therefore  be  ignored  at  the  next  time  step. 
If  the  Reynolds  is  sufficiently  large,  this  error  will  only  affect  the  wake  region 
but  will  not  have  any  influence  at  the  interface  F, .  To  the  contrary,  for  large 
At  and  v,  this  error  does  affect  the  value  of  <j>n+1  on  T*.  The  influence  of 
the  error  on  <j>n+1  may  be  amplified  throughout  the  iteration  process. 

Another  variant  of  the  algorithm  consists  of  replacing  the  explicit  Dirich- 
let  condition 
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=  <f)n  on  rt  in  the  algorithm  (7)-(8) 
by  the  following  semi-implicit  condition 

*Hl  =  r+1  on  Ti. 

In  fact,  this  implies  replacing  the  previously  parallel  algorithm  (Jacobi  like  ) 
by  the  sequential  algorithm  (Gauss- Seidel  like  ). 

When  we  solve  the  pure  diffusion  problem  (i.e.  with  flow  velocity  v  =  0) 
with  v  =  1  and  At  —  1  (respectively  At  =  2)  we  obtain  a  better  convergence 
history  : 

•  the  speed  of  the  new  algorithm  is  linear  and  clearly  faster  than  the 
parallel  algorithm. 

•  the  domain  of  convergence  is  moderately  larger  (see  table  1). 

To  study  experimentally  in  more  details  the  convergence  behavior  of  both 
algorithms  we  assume  that  we  have  a  linear  behavior  of  our  algorithm,  and 
hence  that  the  error  at  the  iteration  n  will  satisfy  the  following  inequality 


Un+l-<t>n\\oo~KnUl-<t>°\\oo 


The  algorithm  converges  if  K  <  1.  An  estimate  for  K  can  be  found  by 
considering  as  in  table  1  the  ratio 


— log 
n 


Un+l  -  <pn\\oo 

||^-^0|ico 


=  -log  K 


which  is  displayed  as  a  function  of  ( vAt )  for  n  =  14  and  different  values  of 

V  —  A  negative  value  of  this  ratio  means  divergence  of  the  algorithm. 
As  expected,  this  ratio  is  positive  for  sufficiently  small  values  of  At  and 
converges  to  zero  as  At  goes  to  zero. 

In  this  table,  we  observe  that  for  V  —  0,  vAt  <  ccq  «  2,  the  algorithm 
converges.  However  the  convergence  is  slow  since  the  minimal  contraction 
constant  Kmin  (for  the  optimal  value  of  i/At )  is  close  to  one  (see  table  2). 
For  V  =  10,  the  algorithm  converges  for  a  much  larger  range  of  values  of  vAt 
and  the  optimal  contraction  constant  is  much  smaller.  This  is  summarized 
on  table  2  where  we  have  displayed  the  best  possible  contraction  constants 
for  each  of  the  coupling  algorithms  and  for  different  values  of  the  Reynolds 

V  =  - 

V  * 
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i/At 

1/10 

1/2 

2 

5 

1000 

Gauss-Seidel 
V  =  0 

0.5 

_ 1 

-0.5 

-0.75 

-0.8 

Jacobi 

V  =  0 

0.1 

-0.09 

-0.25 

-0.4 

-0.41 

Gauss-Seidel 
V  =  10 

0.03 

2.8 

2.6 

2.4 

2.4 

1 

1.15 

1.15 

1.14 

1.14 

Jacobi 

V  =  1000 

0.23 

2.79 

2.8 

2.7 

2.75 

2.8 

Table  1:  Contraction  constant  (in  fact  minus  its  logarithm)  in  function  of 
vAt  for  the  explicit  (Jacobi)  and  semi-explicit  (Gauss-Seidel)  version  of  our 
coupling  algorithm.  We  observe  divergence  for  V  =  0  and  vAt  >  2  and 
convergence  otherwise. 


Jacobi  (parallel) 

Gauss-Seidel  (sequential) 

V  A  min 

V  Kmin 

0  0.85 

0  0.68 

10  0.50 

10  0.11 

IF  044 

Table  2:  Minimal  contraction  constant  versus  the  Reynolds  V  for  both  se¬ 
quential  and  parallel  versions  of  the  algorithm. 
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Figure  2:  Convergence  of  the  Time  Marching  Algorithm:  are 

plotted  versus  the  iteration  count  n  for  uAt  =  10-6,  v  =  0  (Jacobi).  Observe 
the  very  slow  convergence. 
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Figure  3:  Convergence  of  the  Time  Marching  Algorithm: 
plotted  versus  the  iteration  count  n  for  i/A t  =  10_1,  v  =  0  (Jacobi), 


Figure  6:  Convergence  of  the  Time  Marching  Algorithm:  At[[0oJ|co  are 
plotted  versus  the  iteration  count  n  for  vAt  =  10  and  the  flow  velocity  is 
equal  to  1  (Jacobi). 


7  Conclusion 


We  have  analysed  the  convergence  properties  of  a  standard  time  marching  al¬ 
gorithm  for  solving  a  domain  decomposed  advection-diffusion  problem  with 
full  overlapping  and  coupling  by  friction.  We  were  able  to  prove  theoreti¬ 
cally  the  unconditional  stability  and  linear  convergence  of  the  fully  implicit 
algorithm  (§5). 

When  using  the  uncoupled  semi-explicit  algorithm  in  the  general  case, 
numerical  evidence  indicate  that  this  algorithm  is  unstable  for  large  values 
of  At  and  small  overlapping,  and  that  it  becomes  linearly  convergent  when 
At  is  below  a  Reynolds  dependent  threshold  (§7).  This  conditional  stability 
is  not  a  real  issue  for  practical  CFD  problems  because  most  solvers  already 
require  to  use  small  time  steps  inside  each  domain.  Nevertheless,  it  would 
be  nicer  to  derive  an  uncoupled  unconditionally  stable  version  of  the  present 
time  marching  algorithm. 
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Figure  7:  Description  of  the  finite  element  mesh  and  of  the  local  subdomain. 


REPORT  DOCUMENTATION  PAGE 

Form  Approved 

OMB  No.  0704-0188 

. . . . . 

1.  AGENCY  USE  ONLYfieavs  blank)  2.  REPORT  DATE  3.  REPORT  TYPE  AN 

May  1996  Contractor  Repo 

D  DATES  COVERED 

rt 

4.  TITLE  AND  SUBTITLE 

CONVERGENCE  ANALYSIS  OF  DOMAIN  DECOMPOSITION 
ALGORITHMS  WITH  FULL  OVERLAPPING  FOR  THE 
ADVECTION-DIFFUSION  PROBLEMS 

5.  FUNDING  NUMBERS 

C  NAS1-19480 

WU  505-90-52-01 

6.  AUTHOR(S) 

Patrick  Le  Tallec 

Moulay  D.  Tidriri 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Institute  for  Computer  Applications  in  Science  and  Engineering 

Mail  Stop  132C,  NASA  Langley  Research  Center 

Hampton,  VA  23681-0001 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

ICASE  Report  No.  96-37 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

National  Aeronautics  and  Space  Administration 

Langley  Research  Center 

Hampton,  VA  23681-0001 

10.  SPONSORING/MONITORING 

AGENCY  REPORT  NUMBER 

NASA  CR-198339 

ICASE  Report  No.  96-37 

11.  SUPPLEMENTARY  NOTES 

Langley  Technical  Monitor:  Dennis  M.  Bushnell 

Final  Report 

Submitted  to  Mathematics  of  Computation. 

12a.  DISTRIBUTION/AVAILABILITY  STATEMENT 

U  nclassified-U  nlimited 

Subject  Category  64 

12b.  DISTRIBUTION  CODE 

13.  ABSTRACT  (Maximum  200  words) 

The  aim  of  this  paper  is  to  study  the  convergence  properties  of  a  Time  Marching  Alg 
Diffusion  problems  on  two  domains  using  incompatible  discretizations.  The  basic  algoril 
theoretical  or  numerical  results  illustrate  its  convergence  properties. 

orithm  solving  Advection- 
thm  is  first  presented,  and 

14.  SUBJECT  TERMS 

Advection- Diffusion;  Domain  Decomposition;  Time  Marching  Algorithm; 

Spectral  Theory;  Point  Fire  Method 

15.  NUMBER  OF  PAGES 

31 

16.  PRICE  CODE 

A03 

17.  SECURITY  CLASSIFICATION 

OF  REPORT 

Unclassified 

18.  SECURITY  CLASSIFICATION 
OF  THIS  PAGE 

Unclassified 

19.  SECURITY  CLASSIFICATION 
OF  ABSTRACT 

20.  LIMITATION 

OF  ABSTRACT 

NSN  7540-01-280-5500 


Standard  Form  298(Rev.  2-89) 

Prescribed  by  ANSI  Std.  Z39-18 
298-102 


